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Abstract: Fluorescence diffuse optical tomography (fDOT) is an imaging 
modality that provides images of the fluorochrome distribution within 
the object of study. The image reconstruction problem is ill-posed and 
highly underdetermined and, therefore, regularisation techniques need to be 
used. In this paper we use a nonlinear anisotropic diffusion regularisation 
term that incorporates anatomical prior information. We introduce a split 
operator method that reduces the nonlinear inverse problem to two simpler 
problems, allowing fast and efficient solution of the fDOT problem. We 
tested our method using simulated, phantom and ex-vivo mouse data, and 
found that it provides reconstructions with better spatial localisation and 
size of fluorochrome inclusions than using the standard Tikhonov penalty 
term. 
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1. Introduction 

Fluorescence diffuse optical tomography (fDOT) is a relatively new optical imaging modality 
that uses fluorescent markers that accumulate in specific regions, to monitor cellular and subcel- 
lular functional activity. Although it is a technique mostly used for small animal research [1-5], 
it is a promising technique with several medical applications such as detection, diagnosis and 
monitoring of human neoplasms, in particular breast tumours [6-8]. 

In fDOT, a near-infrared (NIR) light source at excitation wavelength is projected onto the 
subject under study at different positions. The excitation light propagates diffusely into the tis- 
sue and some of the photons are absorbed by fluorochromes, which re-emit part of the energy 
at a longer wavelength. Then, fluorescence, and possibly, excitation light intensities are meas- 
ured by detectors placed around the object, for example, opposite the source for transmission 
geometry. This measurement process is repeated for several source positions to obtain a data set 
that is used for the reconstruction of the fluorochrome concentration. In modern fDOT systems, 
detection is performed using a charged-coupled device (CCD) camera, not in direct contact 
with the subject. This leads to large data sets, making the problem largescale and (nominally) 
overdetermined. However, due to the diffusive nature of light propagation in biological tissue 
the image reconstruction is an ill-posed inverse problem [9], hence, in order to obtain a stable 
and meaningful solution one has to use regularisation methods. 

Prior assumptions on the solution, for example its smoothness, can be incorporated in the im- 
age reconstruction [10, 1 1]. Another possibility is to include structural prior information, which 
can be obtained from other imaging methods such as X-ray computed tomography (XCT) [12] 
and Magnetic Resonance Imaging (MRI) [13]. The most common regularisation methods are 
the Tikhonov-type methods, which are L2-norm based and filter high-frequency noise and pro- 
duce smooth images. Tikhonov regularisation with a Laplacian-type operator has been used in 
a few fDOT studies. Davis et al [13] used a Laplacian-type prior that incorporated structural in- 
formation derived from MRI images. Their phantom and simulation studies showed significant 
improvements in the reconstruction of the fluorescence distribution when the structural prior 
was used. Lin et a.l [14] observed that, when both functional and structural Laplacian-type 
prior information is incorporated in the reconstruction, the recovered fluorochrome concentra- 
tion is more accurate. Ale et al. [12] used different Laplacian-type matrices in the regularisation 
functional, and tested the reconstruction method on simulated, phantom and mouse data. Re- 
sults showed that reconstructions were more accurate when their method was used compared 
to the standard Tikhonov method. 
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If the forward model is linearised, by using for example a Born approximation, and the reg- 
ularisation term is L2-norm based, then the reconstruction of the fluorochrome concentration 
is reduced to a linear inverse problem. Nevertheless, in order to preserve edge information in 
the reconstructed image, non-linear regularisation methods, such as total variation (TV) [15], 
can be advantageous. Total variation (TV) regularisation is a LI -norm based method that not 
only penalises highly oscillating solutions, but also preserves edge information in the recon- 
structed image. Image reconstruction using non-linear priors is computationally expensive, and 
convergence can be slow and not always ensured. Freiberger et al. [16] analysed reconstruction 
schemes for fDOT based on the non-linear formulation of the forward model and non-linear pri- 
ors, namely, total variation and level-set type priors. They found, through simulations performed 
using a point source-detection measurement geometry, that the reconstructions performed by 
solving the non-linear problem using a Newton-type method were more accurate than the re- 
constructions obtained using linearised forward models and linear regularisation methods. In a 
previous work, Freiberger et al. [17] introduced an alternating direction minimisation method 
to solve non-linear regularisation methods. This method splits the reconstruction problem into 
a Gauss-Newton step and a TV minimisation step. 

The disadvantage of TV based methods is that in the presence of noise they can lead to the 
staircase effect, i.e., smooth regions appear as piecewise constant regions [18]. Anisotropic dif- 
fusion regularisation with anatomical edge prior has been used in diffusion optical tomography 
(DOT) [19-21] and electrical impedance tomography (EIT) [22] image reconstruction of sim- 
ulated data. This non-linear regularisation technique smooths image noise whilst preserving 
edges, and does not yield staircase effects. 

This paper proposes a split operator method to reconstruct fluorescence images fast and ef- 
ficiently, using anisotropic diffusion regularisation and anatomical prior information obtained 
from another medical imaging modality. The image reconstruction algorithm splits the non- 
linear problem into a linear problem and a nonlinear anisotropic diffusion problem, which are 
straightforward to implement and solve. The performance of the algorithm is assessed using 
simulations, experimental phantom and mouse data. Prior anatomical information is obtained 
from XCT images, which were acquired sequentially with fDOT. Furthermore, the influence 
of using different edge-preserving functions in the anisotropic diffusion prior term is analysed. 
The choice of an edge-preserving function is an important topic, since it defines in which re- 
gions (homogeneous regions/discontinuities) diffusion is isotropic/anisotropic. 

2. Methods 

2. 1 . Formulation of the problem 

In fDOT the continuous -wave forward model is described by a set of coupled diffusion equa- 
tions at the excitation and emission wavelength, X e and Xf, respectively: 

[-V- Jc^A/JV + ^^A/)] U(r,X f ) = U(r,X e )h(r,X f ), r e O (2) 

where *c(r, X) = [3 (ju£(r, X) + /i«(r, A))] -1 is the diffusion coefficient at position r and at wave- 
length A, jl a and jl' s are the absorption and reduced scattering coefficients, respectively. The flu- 
orophore is excited with light at X e emitted by a source q(r,X e ). The emission photon density 
U(r,X e ) over the domain Q is obtained by solving Eq. (1). The fluorescence yield coefficient 
h is related to the quantum yield of the fluorophore 7] and its absorption coefficient at Xf. The 
excitation photon density U(r,Xf) is calculated from Eq. (2). The Robin boundary condition is 
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applied to both equations: 



t/(r,A)+2gir(r,A) X' ' = 0, redQ. (3) 
on 

where £ is a boundary term that incorporates the refractive index mismatch and n is the out- 
ward normal at position r on the boundary <9Q. For weak fluorochrome concentrations the 
contribution of the fluorochrome to the optical parameters can be neglected, which yields the 
first-order Born approximation [23, 24]. Therefore, the measured fluorescence photon density 
yy and excitation photon density y e , due to a source at and a detector at r d , can be written as: 

U(r s ,r d ,he)= / q(r,k e )G(r s ,r d ,k e )dr, 
Jq 

y e = ® e (r s ,r d )U(r s ,r d ,X e ), (4) 



U(r s ,Y d ,X f ) = J h(Y,Xf)U(Y s ,T,ke)G(T,r d ,Xf)dr, 

y f = e f (T s ,r d ,Xf)U(r s ,T d ,X f ), (5) 

where G represents the Green's function solution and 0 are the unknown source and detector 
coupling coefficients. Normalising the measured fluorescence photon density y^ by the meas- 
ured excitation photon density y e reduces the effects of 0 [23, 25]: 

ll = U(r 5 ,T d ,k f ) 
y e U(Y s ,Y d ,X e ) 

= TJ( 1 5 s / h(Y,X f )U(Y s ,YX)G(Y,Y d ,X f )dr, (7) 

where 0/ = Q e is assumed. Note that in CCD based measurement systems, such as used in this 
paper, the previous equation is modified to: 

Jf = P(dQ^E)U(Y s ,Y dl X e ) 
y e P{d£l^Z)U(Y s ,Y d XY 

where the operator P represents the projection from the domain boundary dCl to the camera E. 

The forward problem defined by the set of diffusion Eqs. (1) and (2) is computed numerically 
using the finite element method (FEM) on a tetrahedral mesh based on the geometry being 
considered. However, when using an image based prior, such as the anisotropic diffusion prior, 
it is convenient to map the solution into a regular grid. Therefore, Eq. (6) is discretised into N 
volume elements of size v: 

V 1 N 

J^h{Y U X f )U{Y s ,Y U X e )G{Y U Y d ,Xf)v. (9) 



y e U(Y s ,Y d ,X e ) f = 



For a set of M s source-camera projections, where the detector is a camera with M x x M y 
pixels, the fDOT problem can be written as a linear system: 

y=^=7h, (10) 

where here y^ and y e are vectors of size M s x M x x M y , J is the Jacobian matrix of dimensions 
(M s xM x xM y ) xN, where each row is obtained from Eq. (9) for each source-detector pair, 
and h is a vector of size N. 
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The estimation of the fluorescence yield h from Eq. (10) is often computationally infeasi- 
ble due to the large dimensions of the Jacobian matrix. Therefore, data compression is used 
to reduce the the dimensions of the data for computational efficiency [26], while simultane- 
ously reducing the redundancy of the data. The 2D wavelet transform of a projection image y ; 
captured by a CCD camera can be expanded in the wavelet basis as: 

M x xM y 

Ji= X CiWk (11) 
k 

where c are the wavelet coefficients and <p the wavelet basis. If only the largest My wavelet 
coefficients are kept, which contain most of the relevant information, then the dimensions of 
the vectorised compressed data y are reduced to M s x My. Similarly, the size of the compressed 
Jacobian / is reduced to (M s x M<p) x N. The resulting compressed forward problem is: 

y = /h. (12) 

2.2. Inverse problem methods in image denoising 

Many regularisation schemes in inverse problems originate from image denoising applications. 
One general approach in this application is to convolve an image with a Gaussian kernel, which 
smooths out the noise and edges, and then an anisotropic diffusion filter is applied, which 
recovers the image edges and smooths homogeneous regions in the image [27]. 

2.2.1. Anisotropic diffusion 

The concept of anisotropic diffusion was introduced in image processing by Perona-Malik [28]. 
It is an efficient technique that preserves and enhances edges present in the images while si- 
multaneously removing image noise. Following is a brief analysis of the principle behind this 
method. 

In the image space h E M? the anisotropic diffusion process is given by: 

d -h=V-[g{\Vh\)Vh], (13) 

where \Vh\ = ^/(§) 2 + (f* ) 2 + (§ ) 2 , g(\Vh\) = is the diffusivity, and y is a potential 

function that controls the diffusion process. This term can be written in terms of the gradient 
direction and tangent directions. The unit vector in the direction of the gradient is 7] = and 
the two unit vectors in the tangent direction, i.e., orthogonal to the gradient, are ^\ and <^2- Using 
the divergence product rule of a scalar function / and vector v: V- (/v) = / (V- v) + V/- v, the 
anisotropic diffusion operator can be written as: 



W(lv*l). VA 



v/(|v/i|)v '(^l) +V(v/(|V,i|)) ^- (14) 



After some manipulation we have the directional anisotropic diffusion: 



|VA| 



+H^) +V"(I VA|)A,„ (15) 



where h^^, and are the second derivatives of h in the £i, §2 and Tf directions, 

respectively. This formulation can help us understand the effect of the anisotropic diffusion 
regularisation in the normal and tangential directions. In order to preserve edges the regularisa- 
tion should locally penalise magnitudes of the directional derivative across an edge much less 
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than along an edge. Thus, the function \f/ should be constructed such that the diffusion along 
the normal direction r\ is small compared to the tangential directions and <^2- 

The isotropic diffusion (Tikhonov) regularisation, i.e., when i/a(|V/z|) = ^|V/i| 2 , does not 
preserve edges since it acts in all directions: 

V(|VA|) 



-Vh 



: + A «i«2 + A n»i- (16) 



2.2.2. Discretisation 



The discrete formulation of the anisotropic diffusion in Eq. (13) using an explicit scheme, at 
the k th iteration, is given by [29] : 

UM uk o k + o k 



h^ l -h k _ 



sP-.I^fe-*?)- (17) 



where gi is the diffusivity of pixel i, K (i) are the neighbours of pixel i, and At in the time step 
and d is the discretisation interval. For m-dimensions, the previous equation can be written in a 
matrix-vector notation as [29] : 

LA = ^^ ( ^ ( i8) 

Ar 1=1 

h M = |/ + A* X-2?z(A*)^ A*, (19) 

where the matrix <S?i(h k ) = [Uj{h k )\ corresponds to the diffusivity derivatives along the /th 
coordinate axis: 



'41 jeK(0 



0 ^/i 1 ^ 



Note that only two neighbours of pixel i are considered for each direction. The drawback of the 
explicit scheme is that it is only conditionally stable. This means that the step size needs to be 
small for the solution to be stable, and therefore, a large number of iterations are required. For 
example, in 2D At < 1/4. The allowed step size becomes smaller for higher dimensions. 

Alternatively, we use the semi-implicit additive operator splitting (AOS) method, which as- 
sumes the following expression [29] : 

1 m / \ -1 

h M = _ £ h- m At^f l (h k )) h k , (20) 

This scheme is unconditionally stable. Thus, At can be relatively large without causing nu- 
merical instability. Since it is an additive method one can treat each direction separately. 
Therefore, the AOS scheme involves solving m linear systems Dh k+l = h k , where the ma- 
trix D = (/ — At mJ??i(h k )) is tridiagonal. Note that for each direction there is one tridiagonal 
matrix. An efficient way of solving this is to use the Thomas algorithm, which is a Gaussian 
elimination algorithm for tridiagonal systems [29] . This method avoids the need of directly de- 
termining the inverse of large matrices, since the problem is reduced to computing a series of 
multiplications, divisions, additions and subtractions. 
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Table 1 . Edge-Preserving Functions 



Perona- Malik (Cauchy or Lorentzian) [30] 


y(|Vft|) = £io g ( 


2 




Perona- Malik 2 (Welsh) [30] 


y(|VA|) = £ 
g(|Vfc|)=exp 


1- 

(- 


exp | 

( ™ ; 






2 )] 




Total variation approximation [21] 


V (|Vfc|) = r v 
*(|v*l) = ^ 




-T 




r 

/i| 2 +r 2 


Huber [21,30] 


yf(\Vh\) = < 
*(|VA|) = | 


f r|VA|-2 
I ¥ 

r 


\Vh\ < T 
\Vh\ > T 
\Vh\ < T 




Tukey [30] 


¥(|V/i|)=< 
«(|VA|) = | 


[ 


r 2 

6 

T 2 
6 

1- 

) 


1- 


1- 

r 


2 

1 

r 


)1 

^| 


3" 

\Vh\ < T 
\Vh\ > T 

< T 
> T 


Exceedance 


y(\Vh\) = - 
g(|VA|) = l 


¥f-(l-P(\Vh\ <X)) 
-P(\Vh\<X)=P(\Vh\>X) 



2.2.3. Edge-preserving functions 

As mentioned previously, the diffusion regularisation should be isotropic in homogeneous re- 
gions and in the presence of an edge the diffusion regularisation should be significantly smaller 
in the normal direction than in the tangential direction. Table 1 shows a few functions with 
these edge preserving properties [21, 30]. In the isotropic diffusion case we have the first-order 
Tikhonov regularisation, where \j/(\Vh\) = ^|V/z| 2 andg(|V/i|) = 1. 

The parameter T is the threshold. If this parameter is too large it leads to an oversmoothing 
of the image. Whereas, if T is too small, smoothing is not applied and the resulting image is 
similar to the initial one. The parameter T can be selected using a method based on normalised 
cumulative histogram (NCH) of the gradient [28]. This method is commonly used in edge de- 
tection problems. The NCH indicates the probability P of a gradient taking on a value less 
than or equal to the value X that the bin represents. It increases monotonically and the smooth- 
ness/sharpness of the curve indicates how smooth/sharp the edges are. The threshold can be 
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calculated from the NCH by setting the threshold at, for example, 90 percent. It can be kept 
fixed or it could be updated at every iteration. Note that for each edge-preserving function the 
parameter T is scaled, so that all functions preserve edges at the same threshold. 

The exceedance function can be thought of as the probability of an edge of interest being 
present and it can be calculated using the NCH of the image gradient. The NCH gives the 
probability P(\Vh\<X). Consider X = T, then P(\ Vh\ < X) ~ 0 and g ~ 1, on the other hand 
NCH tends to 1 as the bin values increase, and therefore, g tends to 0. The exceedance function 
has a similar behaviour than the other edge-preserving functions, but does not require a direct 
estimation of T. 

2.3. fDOT inversion using anisotropic diffusion regularisation 

We now define the inverse problem of fDOT as minimising the following objective function: 

1 ,, ~ ,,2 

Minimise E(h) = - \\f-Jh\\ + a*F(h), (21) 

where a is the regularisation parameter, and *F(h) is the prior functional given by: 

*F(h) = / v(|Vh|)dO. (22) 

The minimum of E(h) can be calculated using the gradient: 

V£(h) = f (Jh - y) + a JSf (h)h, (23) 
where JSf is a linear operator given by: 

JS?(h) = -V.b(|VA|)V]. (24) 
2.3.1. Incorporation of an anatomical prior 

If prior information about the image structure is available, it can be used in the prior term to 
weight the diffusion function. Consider x re f to be an anatomical image, which has been filtered 
previously by a Gaussian or anisotropic diffusion filter. Then, in order to obtain the image 
weighting factor W(|Vx r£? /|) one can apply one of the diffusion functions in Table 1 to the 
anatomical image, which will return an image where homogeneous regions have a large weight 
and close to edges the weight approaches zero. The prior term takes the form: 

V w (h) = W(|V^/|)v(|Vh|)dQ, (25) 

which leads to 

Jzfw(h) = V- [W(\Vx ref \)g(\Vh\)V] . (26) 

Instead of keeping the threshold parameter T (see §2.2.3 ) constant over the entire image 
we use a spatially variant T, which is updated at each iteration. The normalised image gradient 
V/i* and normalised structural weighting factor W^(|x re /|) can be used to weight the threshold 
value found from the NCH. Note that both quantities are normalised so that they are equal 
to 1 in flat regions and tend to 0 near the edges. The spatially variant threshold is calculated 
as: T(x,y,z) = T + T [g n (x,y,z) — 1/2)], where T is the threshold found using the NCH and 
g n (x,y,z) = (yh„ + W n (\x re f\)} 1 2. Therefore, in the presence on an edge in either the optical 
image or the anatomical image, the parameter T is small and the edges are preserved, since the 
diffusion process is stopped. In homogeneous regions T is large, which results in a stronger 
smoothing effect. The variant T is particularly useful for restoring features in the image with 
different contrasts. 
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2.4. Image reconstruction algorithm 

The Gauss-Newton approach can be used to obtain the solution iteratively [31]: 

h* +1 = h* + t (j T J + a*F"(h*)) _1 f ( (y - M^j - aY(h k )^j , (27) 

where T is the step length. In the lagged-diffusivity Gauss-Newton method the Hessian of 
the prior is calculated from a linear diffusion operator: ^"(h) = j£f (h). However, the Gauss- 
Newton method Eq. (27) requires a computationally expensive matrix inversion. In order to 
overcome this issue and aiming at having fast runtime performance, a split operator algorithm 
is introduced in this paper. 

We want to find an update step for Eq. (21). We note that Jzf (h) is a symmetric positive 
definite regularisation matrix, and we assume both a decomposition and inverse to exist such 
that: 

L T L = jSf (h), r = [J^(h)]- 1 = L-\L- l ) T . (28) 

Since Jzf (h) is a differential operator, then F is a smoothing operator. Now consider the trans- 
formation: 

J = JL-\ h = Lh, (29) 

such that Eq. (21) becomes: 

Minimise E(h) = ^ ||y-/h|| 2 + a \\h\\ 2 . (30) 

An iterative descent method (the Landweber method) [32] for the previous regularised least 
squares problem (Eq. (30)) is: 

h M =h k -F(y + Jh k ) (31) 
h M =h k + Tf (y - /V) . (32) 

We can interpret Eq. (31) as a gradient descent step for the likelihood (first term is Eq. (21)), 
followed by a smoothing step to regularise the update: 

h w / 2 = h^+/(y-^) (33) 

^^+1)1^+1^+1/2 (34) 

In order to achieve a faster convergence we may use a Levenberg-Marquardt (LM) step instead 
of a descent step: 

h k+l/2 =h k + JT B ^ _ jtf>j ^ (35) 

where B is a symmetric positive definite matrix. Here we define: 

B = (JJ 7 + (36) 

which is cheap to invert because we compressed the row space of the Jacobian using wavelets. 
The parameter A in the LM method is the damping factor, which changes at each iteration 

Full inversion of the smoothing step is expensive, so we replace the regularisation step by 
the semi-implicit AOS method (Eq. (20)), which can be solved efficiently using the Thomas 
algorithm. We may now relax the requirement that an inverse of J2?(h) exists, and we do not 
need to consider the decomposition used in Eq. (28). 
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A simpler way to infer this split operator method is by considering Eq. (23) as the update Ah^ 
of the gradient descent method. Since Ah^ consists of two terms it can be split into a two-step 
iteration. 

This reconstructions method involves solving two systems of equations, where the first step 
is simply the LM method [32] and the second step is the semi-implicit form of the anisotropic 
diffusion method, where At is the time step, which plays the role of a in Eq. (21) [33]. There- 
fore, the first step reconstructs the fluorescence distribution image and the second step applies 
the prior, by deblurring/denoising the reconstructed image. 

The parameter A in the LM method serves a role similar to the step size T, since it controls 
the magnitude and direction of the iteration update Ah k . The LM method switches between the 
Gauss-Newton method and the gradient descent method. When A —> 0 the method is basically 
the same as Gauss-Newton. On the other hand when A —> <*> the method becomes a gradient 
descent approach. Since the LM method is a special case of Tikhonov regularisation [34, 35] 
and, furthermore, h° = {0, ..,0}, then A 0 can be calculated by the L-curve method A/ c [36, 37]. 

For the following iterations, if the least squares norm 

y_y h £+i/2 < ||y-/V|| thenA* +1 
is decreased by a certain percentage of the current value, px , since the Gauss-Newton method 



converges quickly near a minimum. If 



> 1 1 y — Jh k 1 1 then A kJr 1 is greater than 



X k , since the gradient descent method is more efficient when the procedure is far from the 
minimum. 

The image reconstruction is summarised in the following algorithm: 



Algorithm 1 Two-Step Image Reconstruction Algorithm 



a 0 = A/ c trace (JJ T ) 
a = 0 

while ||y-./h*|| 2 <eW k<N k do 

h *+l/2 =h k + jT (JJT + A/ ) - 1 (y _ ytf) 

for i= 1 , , n do 
for / = 1 , , m do 

v= (/ + mA/^(h w /2)) _1 h Hi/2 

u = u-\- V 
end for 

m 

end for 



e = tolerance value 



AOS -Thomas algorithm 



if y-yh* +1 / 2 

X k + l =X k -X k * Px 
else 

end if 
end while 



< Hy-Jh^ll 2 then 



From our experience, for the prior step, it is preferable to perform a small number of iterations 
n and use a relatively large At, than just a single iteration and a very large At. 
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2.5. Evaluation 

The performance of our image reconstruction algorithm is evaluated through simulations, ex- 
perimental phantom measurements and ex — vivo mouse measurements. 

The image reconstruction algorithm implementation is based on the software package Time- 
resolved Optical Absorption and Scattering Tomography (TOAST) developed at University 
College London (UCL). TOAST is a finite element method (FEM) based software and therefore 
requires a suitable computational mesh. Meshes were generated from the XCT images of the 
phantom and mouse. These images were contaminated by noise and were smoothed using a 
simple Gaussian filter. The images were segmented by thresholding, and the segmented images 
were used to generate a tetrahedral mesh using the Matlab package Iso2mesh [38]. 

The simulation procedure and experimental setup are described below. A tolerance value £ = 
le— 4 and maximum iteration number A/& = 30 were used as the stopping criteria in Algorithm 1 . 

2.5.1. Simulations 

The Digimouse atlas [39] was used to generate a mouse mesh with 76 973 nodes, 450 829 
elements and dimensions 32.5 mm x 19.5 mm x 88 mm. A fluorescent target was used to sim- 
ulate a tumour in the liver. The target is a sphere of radius 1.75 mm and contrast h = 1 mm -1 
(h = 0 mm -1 outside). We used a simplified two-tissue model, where different optical proper- 
ties were assigned to the liver (fi a = 0.035 mm -1 and fl' s = 0.68 mm -1 ) and other tissue (\l a = 
0.01 mm -1 and fl f s = 0.8 mm -1 ). The optical properties were calculated using the parameters 
provided by Alexandrakis et al [40] at a wavelength of 670 nm. The source and detector (CCD 
camera), which was placed opposite the source, were considered to rotate around the mouse. 
Projections were calculated for 16 evenly-spaced positions over the full 360°range. Data con- 
sisted of fluorescence and excitation projections with 1% additive Gaussian random noise. For 
the reconstruction a total of 128 wavelet coefficients were used. 

The peak signal-to-noise ratio (PSNR) in decibels (dB) is used to evaluate the quality of the 
reconstructed images. The PSNR between the true image h true of dimensions X xY xZ and the 
reconstructed image h recon is defined as follows: 

™=101og 10 ^pl, (37) 



where MSE is the mean squared error: 



MSB = g^jkjk^g h ™°J\ 
XxYxZ 

A higher PSNR indicates a greater fidelity to the original image. 
2.5.2. Experimental setup 

2.5.2.1. fDOT-XCT system The fDOT-XCT rotating system was developed and assembled 
at Universidad Carlos III de Madrid (Spain) [41]. The fDOT and CT components were mounted 
on a rotating gantry, whose rotation movement was controlled by a Newport RV350PP mo- 
torised rotation stage (Newport Co, Irvine, CA, USA). The object to be imaged was placed in 
the horizontal position on a bed specifically designed for the experiment purposes, which was 
mounted on a motorised linear translation stage (LTM80-300-HsM, OWIS GmbH, Staufen, 
DE) to accurately position the object in the field of view (FOV) of the system, and was gently 
compressed between two methacrylate transparent plates to produce planar surfaces. 

Note that in our studies the gantry did not rotate during the experiment, and only the source 
position changes. The CCD camera (ORCA II, Hammamatsu, Japan) and objective lenses 
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(Nikon) was placed parallel to the anti-reflective plates. The laser diode (Monocrom, Barcelona, 
Spain) with an emission wavelength of 675 ± 5 nm was placed opposite the CCD camera. The 
output power was controlled via pulse width modulation (PWM) using transistor - transistor 
logic (TTL) levels. The laser light was focused onto the sample at predefined points via two 
mirrors moved by galvanometers (Scancube 7, ScanLab AG, Puchheim, Germany). For ev- 
ery source position, the laser power that gives an optimal number of counts recorded by the 
CCD camera was obtained using an automated algorithm. Fluorescence measurements were 
recorded by placing a 10 nm bandwidth filter centred at 700 nm in front of the objective of 
the CCD camera, whereas for excitation measurements a 10 nm bandwidth filter centred at 675 
nm was used. A motorised filter wheel (Luxiflux V2, Cyberstar, Echirolles, France) was used 
to quickly switch between filters. All the acquisition processes were controlled via a software 
written in C++ and IDL code. 

The micro-CT consisted of a modified Hamamatsu L963 1 microfocus X-ray source (Hama- 
matsu Photonics, Japan), in which the X-ray control unit was removed from its housing in order 
to reduce size and weight of the device attached to the rotating gantry, and a flat-panel X-ray 
detector (Hamamatsu C7940DK-02, Hamamatsu Photonics, Japan). The X-ray source has an 
output power of 50 W, a maximum peak energy of 1 10 keV and a maximum anode current of 
0.8 mA. The detector pixel size was 0.05 mm and up to 4x4 pixels could be combined through 
pixel binning, in order to speed image acquisition (at the expense of image resolution). The 
image integration time was 125 ms when 4x4 binning was used. 

The detector was mounted on a motorised linear stage (LTM80-100-HSM, OWIS GmbH, 
Staufen, DE), which allows the X-ray detectors to be moved radially to set the desired FOV 
size and magnification factor (that determines the resolution of the acquired images). Images 
were acquired with a fixed integration time of 125 ms, for a total of 360 projections over an 
angular range of 360°. The resulting images had a voxel size of (0.144 mm) 3 . 

2.5.2.2. Phantom A slab phantom made of a polyester resin was built with ji a = 0.01 mm - 1 
and \±' s = 0.8 mm -1 [42]. The optical properties of the phantom were considered to be the 
same at both measurement wavelengths. The phantom dimensions were 50x50x10 mm. A 
rectangular mesh with slighter larger dimensions than the phantom was generated (15x15 x 10 
mm), which contained 65 596 nodes and 303 750 elements. A capillary of 1 mm diameter was 
inserted into the phantom, close and parallel to the imaging surface, and its tip was filled with 
2 jiL of Alexa Fluor 680 with fluorochrome concentration 20 [i Molar (jiM). The capillary was 
located approximately 1.8 mm below the imaging surface and at the centre (x = 25 mm) of the 
phantom. The extremities of the capillary tube were at y = 0 mm and y = 31 mm. A total of 
42 projections were obtained (the CCD camera, placed opposite the source, was static and only 
the source position varied). Sources were equally spaced within a square region of 10 x 10 mm. 
Images were reconstructed with 128 wavelet coefficients. 

2.5.2.3. Mouse Similar to the phantom experiment, a capillary with 2 jiL Alexa Fluor 680 
with fluorochrome concentration 20 /iM was inserted into the esophagus of an euthanised adult 
nude mouse. A total of 56 equally spaced sources scanned a 10x10 mm region of the upper 
torso of the animal. A finite element mesh was generated, from the XCT images, with 115 764 
nodes and 608 319 elements. For simplicity, optical properties were assumed to be homoge- 
neous (jl a = 0.01 mm -1 and /A f s = 0.8 mm -1 ) and the same for both excitation and emission 
wavelengths. A total of 128 wavelet coefficients was used in the data in the image reconstruc- 
tion. 
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Fig. 1. Reconstructions of h (mm -1 ) from simulated data (image dimensions 32.5 mm x 
88 mm) and PSNR values (dB): (a) the target, (b) reconstruction with Tikhonov prior, (c) 
reconstruction with anatomical prior and total variation, (d) reconstruction with anatomi- 
cal prior and Perona-Malik function, (e) reconstruction with anatomical prior and Perona- 
Malik 2 function, (f) reconstruction with anatomical prior and Huber function, (g) recon- 
struction with anatomical prior and Tukey function, and (h) reconstruction with anatomical 
prior and exceedance function. 

3. Results 

3.1. Simulations 

The simulated fluorescent inclusion in the mouse liver and reconstructed images of the fluo- 
rescence yield coefficient (in mm -1 ), overlay ed on the anatomical atlas, are shown in Fig. 1. 
For comparison with the images reconstructed using the proposed method, Fig. 1(b) shows 
the image reconstructed using a simple Tikhonov regularisation. Figures l(c)-l(h) show the 
reconstructed images using the anatomical prior and different edge-preserving functions, with 
spatially variant threshold T. The PSNR value for each reconstruction is also displayed. Im- 
ages reconstructed using the Perona-Malik 2, Tukey and exceedance functions have the highest 
PSNR values (in decreasing order). 

Figure 2 shows the weighted Perona-Malik edge prior at the first and last iteration of the 
reconstruction. The prior at the first iteration only shows the anatomical edges. As the number 
of iterations increase, the edges of the liver weight the influence of the edges in the reconstructed 
image. Therefore, diffusion is blocked across boundaries where the weighted prior has a small 
value, thus enhancing the edges of the reconstructed image, as well as smoothing the remaining 
regions where the weighted prior has a larger value. 

Figure 3 shows the normalised profiles across the fluorescent target, for the longitudinal and 
lateral directions, for the images reconstructed using Tikhonov regularisation and our method 
with the edge-preserving function that returned the highest PSNR value. 

3.2. Phantom 

The phantom consisted of a resin slab with a small capillary tube filled with a fluorescent dye, 
which was located close to the imaging surface. The XCT images were used to extract the cap- 
illary structure, i.e., the structural prior W n (|jc rg /|). Figure 4(a) shows the fluorochrome concen- 
tration (in /iM) reconstructed using Tikhonov regularisation. The reconstructed concentrations 
using the structural prior and the different potential functions are shown in Figs. 4(b)-4(e). Fig- 
ure 4(b) shows the reconstruction with total variation, Fig. 4(c) shows the reconstruction with 
Perona-Malik function, Fig. 4(d) shows the reconstruction with Perona-Malik 2 function and 
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(a) (b) 



Fig. 2. Weighted Perona-Malik prior (a) at the first iteration, where only the anatomical 
edges are visible, and (b) at the last iteration, where both anatomical and fluorescence 
target edges can be seen. 




x(mm) " J " y(mm) 



(a) (b) 



Fig. 3. Normalised profile plots across the fluorescent target: (a) lateral direction and (b) 
longitudinal direction. 




(a) (b) (c) (d) (e) 



Fig. 4. Fluorochrome concentrations (fiM) reconstructed from phantom data using : (a) 
Tikhonov prior, (b) structural prior and total variation, (c) structural prior and Perona- 
Malik function, (d) structural prior and Perona-Malik 2 function, and (e) structural prior 
and exceedance function. Images have dimensions 150 mm x 150 mm. 

Fig. 4(e) shows the reconstruction with the exceedance function. Similar reconstruction were 
obtained using the Huber and Tukey methods, and therefore are not shown. 

3.3. Mouse 

The XCT images of the mouse with a capillary inserted in the esophagus were used to obtain the 
anatomical prior, where the main visible structures were the bones and the capillary. Figure 5 
shows the reconstructed images of the fluorochrome in the capillary, overlayed on the XCT 
image. Figure 5(a) shows the recovered fluorochrome concentration (in /iM) obtained using 
the Tikhonov method and Fig. 5(b) shows the image reconstructed using the Perona-Malik 2 
method weighted by the anatomical prior. 
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Fig. 5. Fluorochrome concentrations (fiM) reconstructed from mouse data using : (a) 
Tikhonov prior, (b) structural prior and Perona-Malik 2. Images have dimensions 32 mm 
x 32 mm. 

The computational time of a single iteration was approximately 1.7 s on a 2.66 GHz Intel 
dual core processor. The anisotropic diffusion step was computed using the AOS -Thomas algo- 
rithm in 0.92 s. It took approximately 70 s to run the algorithm until the stopping criterion was 
reached. However, the reconstructions obtained after the first 5 iterations were already quite 
good, which is achieved in less than 9 s. 

4. Discussion and conclusion 

In this work we introduced an operator splitting method for solving the fDOT inverse prob- 
lem using non-linear anisotropic diffusion regularisation. A Levenberg-Marquardt minimisa- 
tion step alternates with an anisotropic diffusion filtering step, which can be thought of as an 
image reconstruction step followed by an image denoising/deblurring step. We presented a fast 
and efficient algorithm for the anisotropic diffusion step, based on the semi-implicit additive 
operator splitting method. The performance of our reconstruction method was tested on simu- 
lated, experimental phantom data and ex-vivo mouse data. We also analysed the performance 
of a few different edge-preserving functions that had been proposed in the literature. It is im- 
portant to choose an adequate edge-preserving function g and threshold value T to be used in 
the anisotropic diffusion prior, since it determines which gradients are true outliers, and hence, 
if smoothing should be inhibited in order to preserve edges. We have introduced a T that varies 
according to the gradient of the reconstructed optical image and anatomical weighting factor, 
and is updated at each iteration. 

Simulations showed that our proposed splitting method with anisotropic diffusion regu- 
larisation provides better reconstructions than the standard Tikhonov method. The different 
edge-preserving functions have very similar effects on the reconstructions, since they all are 
monotonically decreasing functions, except the exceedance functional, which penalise gradi- 
ents smaller than the threshold value, and for larger values these functions are almost flat and 
tend to zero, meaning that these gradients are preserved in the images. However, we observed 
that edge-preserving functions that are zero (or very closer to zero) for the gradients considered 
to be outliers (Tukey and Perona-Malik 2 functions) provide better localised images, which are 
also qualitatively more accurate, than other functions that, even though in a small amount, still 
penalise the outliers (Huber and total variation functions). The exceedance functional has the 
advantage that it does not depend on the parameter T, since it is based on the probability of a 
gradient being an edge. 

In the phantom study, even though the capillary was placed quite close to the imaging surface 
and the reconstruction obtained using the Tikhonov prior is already quite good, it is clear that 
our method with the anisotropic diffusion prior provides better images by removing the noise 
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in the homogeneous regions and enhancing the capillary edges. In the mouse study, the effect 
of the anatomical prior on the reconstruction is quite evident. The anatomical prior blocks dif- 
fusion across the edges, hence the reconstructed fluorochrome concentration appears confined 
within the capillary boundaries. 

Our results show good estimates of the dimension and location of the fluorophore distribu- 
tion, which was the main goal of this work. This information is essential in applications such 
as cancer treatment using radiotherapy, so that radiation is directed towards the tumour while 
healthy surrounding tissue receives a minimum amount of radiation possible. It is important to 
note that the fluorochrome is not visible in the reference anatomical image. 

An important feature of our method is that it reflects both structural information present in the 
reference image and also that present in the fluorescent image. Reconstruction do not constrain 
the fluorophore to be confined by the structures in the reference image, but allow for it to be 
distributed across tissue boundaries. 

The computation of the Jacobian can be time consuming depending on the number of 
measurements, number of wavelet coefficients used and number of mesh nodes. For exam- 
ple, it took 165 s per source to built the Jacobian for our simulation study. This calculation can 
be twice as fast if half the wavelet coefficients are used for the data compression. For a mesh 
with 20% of the nodes of the mesh used in our simulations and 128 wavelets, it takes about 
17 s. Alternatively, the Jacobian matrix could be calculated using a matrix-free algorithm [11] 
or by parallelisation of the Jacobian calculation. 

The underestimation of h observed in the simulations needs to be further investigated. A 
possible explanation for these results is the partial volume effect, which is known to lead to 
an underestimation of the fluorescence yield. Furthermore, the simulated target contrast in the 
nodal basis is lower than in the voxel basis. 

In summary, our operator splitting method is straightforward to implement and is compu- 
tationally fast. The images reconstructed using the anisotropic diffusion prior with anatomical 
prior show superior quality than those reconstructed by solving a simpler linear problem using 
Tikhonov regularisation. The weighted anisotropic diffusion prior combines the advantage of 
preserving/enhancing the edges of the reconstructed optical image while reducing the noise and 
stopping diffusion through boundaries defined by the anatomical prior. However, further inves- 
tigation is required to evaluate the adequacy of our method for fDOT image reconstruction of 
in-vivo animal data. 
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